% Michael Sieber, Philipp Rusch

function[alpha, stats] = PolynomInterpolation(lbound, rbound, fun, x, s, kmax)
	counter = 1; % iteration counter    
	
	alpha0 = 0; % left bound
	alpha1 = 1.75; % right bound
	alpha2 = 1; % approximated optimum
	
	epsilon = 10^-4;
	
	%iteration loop
	do
		% equation matrix
		A = [alpha0^3, alpha0^2, alpha0^1, 1; 
			 alpha1^3, alpha1^2, alpha1^1, 1; 
			 3*alpha0^2, 2*alpha0^1, 1, 0;
			 3*alpha1^2, 2*alpha1^1, 1, 0];
		
		% calculate phi
		[f,g] = feval(fun,x + alpha0*s);	
		[k,l] = feval(fun,x + alpha1*s);
		% approximation vecotr 
		b = [f, k, g'*s, l'*s]';
		% calculate constants
		y = A\b;
		
		
		% solve equation
		p = quadsolver((2*y(2))/(3*y(1)), y(3)/(3*y(1)));
		% get new alpha values
		if(p(1) > alpha0 && p(1) < alpha1)
			alpha2 = p(1);
		elseif(p(2) > alpha0 && p(2) < alpha1)
			alpha2 = p(2);
		endif
		
		% set new alpha values
		%if(gradient(alpha2) > 0)  gradient(feval(fun,x+alpha2*s))
		[m,n] = feval(fun,x+alpha2*s);
		
		% steigung in suchrichtung s
		steigung = n'*s;
		if(steigung > 0)
			alpha1 = alpha2;
		else
			alpha0 = alpha2;
		endif

		% some stats
		alpha(counter) = alpha2;
		stats(counter) = counter;
		
		% calculate new phi
		[fnew,gnew] = feval(fun,x + alpha2*s);
		counter = counter + 1;
		gnew'*s
	until (counter > kmax || abs(gnew'*s) <= epsilon)
	%x+alpha2*s 
endfunction % end of polynomial interpolation